In vivo transomic analyses of glucose-responsive metabolism in skeletal muscle reveal core differences between the healthy and obese states

Metabolic regulation in skeletal muscle is essential for blood glucose homeostasis. Obesity causes insulin resistance in skeletal muscle, leading to hyperglycemia and type 2 diabetes. In this study, we performed multiomic analysis of the skeletal muscle of wild-type (WT) and leptin-deficient obese (ob/ob) mice, and constructed regulatory transomic networks for metabolism after oral glucose administration. Our network revealed that metabolic regulation by glucose-responsive metabolites had a major effect on WT mice, especially carbohydrate metabolic pathways. By contrast, in ob/ob mice, much of the metabolic regulation by glucose-responsive metabolites was lost and metabolic regulation by glucose-responsive genes was largely increased, especially in carbohydrate and lipid metabolic pathways. We present some characteristic metabolic regulatory pathways found in central carbon, branched amino acids, and ketone body metabolism. Our transomic analysis will provide insights into how skeletal muscle responds to changes in blood glucose and how it fails to respond in obesity.


Results
Overview of the study. Metabolic reactions, which are defined as chemical reactions of metabolism, are regulated by an integrated network of metabolites as allosteric regulators, substrates, and products; metabolic enzymes; transcription factors; and signaling molecules. To elucidate the regulatory network controlling glucose-responsive metabolic reactions in skeletal muscle, we constructed a regulatory transomic network by integrating metabolic reactions with metabolites, gene expression of metabolic enzymes, and transcription factors, using skeletal muscle excised from C57BL/6J WT mice or ob/ob mice at different time points after glucose administration (Fig. S1). The transomic network of the skeletal muscle was constructed according to our previous study of the liver 12 .
Glucose was administered orally to 16 h-fasted WT and ob/ob mice, and the gastrocnemius muscle and blood were collected at 0, 20, 60, 120, and 240 min after glucose administration (Fig. 1A). The ob/ob mice showed elevated levels of blood glucose and insulin compared to WT mice throughout the study (Fig. S2A). The blood and skeletal muscle data in the fasting state were obtained from our previous studies 12,13 . The skeletal muscle data after oral glucose administration were newly obtained in this study (Fig. S2B).
Using the skeletal muscle data during oral glucose administration, we defined the features of glucose-responsive molecules according to our previous study 12 . Molecules that showed statistically significant changes (absolute log 2 fold change ≥ 0.585 [2 0.585 = 1.5] and a false discovery rate [FDR]-adjusted p value [q value] ≤ 0.1) at any time point compared to the fasting state after glucose administration were defined as glucose-responsive (Fig. 1B). We also calculated time constants (T 1/2 ) to study the temporal patterns of glucose-responsive molecules (Fig. 1C). T 1/2 was defined as the amount of time needed for the response to reach half of the minimum (decreasing molecules) or maximum (increasing molecules) amplitude. According to the blood insulin concentration, which peaked at about 20 min and decreased to basal level at about 60 min (Fig. S2A), rapid responses were defined as those with T 1/2 values less than 20 min, and slow responses were defined as those with values longer than 60 min.
Glucose-responsive molecules were integrated across the omic layers, and the regulatory transomic network was constructed in WT and ob/ob mice (Fig. 1A). The transomic networks contained layers of insulin signaling molecules (Insulin signal), transcription factors (TF), gene expression and phosphorylation of metabolic enzymes (Enzyme), metabolic reactions (Reaction), and metabolites (Metabolite), and the layers were connected when regulations could be speculated. By comparing the regulatory transomic networks between WT and ob/ob mice, we comprehensively evaluated how obesity affects the responses to glucose in skeletal muscle.

Metabolomics.
We first performed metabolomics analysis using capillary electrophoresis-mass spectrometry (CE-MS), liquid chromatography (LC) -MS, and enzyme assays. A total of 104 water-soluble and ionic metabolites including glucose, amino acids, and nucleic acids were measured by CE-MS. Statistical tests were www.nature.com/scientificreports/ performed to identify the glucose-responsive metabolites in WT and ob/ob mice ( Fig. 2A,B; Data File S1). To define an increase or decrease in time courses with changes in both directions at different time points, the direction of change compared to time 0 at the earliest time point that showed a significant change was used.
Metabolites that showed statistically significant increases or decreases in WT or ob/ob mice are shown in Fig. 2A.
The responses were categorized into three groups (rapid, intermediate, or slow) according to their T 1/2 values (Fig. 2C). Four metabolites (4% of the total quantified metabolites) were significantly increased only in ob/ob mice, and none were increased in WT mice (Fig. 2B). Metabolites that were increased only in ob/ob mice included fructose 6-phosphate (F6P), tartrate, lactate, and isoamylamine (Fig. 2D). Twenty-one metabolites (20%) were significantly decreased only in WT mice, and four metabolites (4%) were significantly decreased only in ob/ob mice (Fig. 2B). It is noteworthy that no common metabolites were increased or decreased in WT and ob/ob mice. Metabolites decreased in WT mice included those that play a role in the tricarboxylic acid (TCA) cycle, such as citrate, cis-aconitate, succinate, and malate (Fig. 2D). The ketone body 3-hydroxybutylate (3-OH-butylate) was also decreased in WT mice. Metabolites that were decreased in ob/ob mice included valine, aspartic acid, choline, and 3-methylhistidine. Most of the decreased metabolites showed rapid responses in both WT and ob/ob mice (Fig. 2C). Hierarchical clustering analysis of the metabolites is shown in Figure S3. LC-MS did not detect significant responses of 14 lipids after oral glucose administration (Data File S2).
Our metabolomic analysis revealed that the number of glucose-responsive metabolites specific to WT mice (21: 0 increased + 21 decreased) was larger than that specific to ob/ob mice (8: 4 increased + 4 decreased), and no responses were common to both mice. These results indicate that there is a substantial difference in the mechanism of glucose metabolism in skeletal muscle between WT and ob/ob mice.
Next, we compared the metabolomic changes in the skeletal muscle and blood. The amount of metabolites was regulated not only within each organ but in the blood circulatory system 14 . For each metabolite that was measurable in both skeletal muscle and blood (61 metabolites), we calculated the correlation between the time course of Figure 1. Pipeline of the construction of the glucose-responsive transomic network. (A) We measured the time courses of multiomic data from the skeletal muscles of WT and ob/ob mice following oral glucose administration and identified the molecules that were changed by oral glucose administration, which we defined as glucose-responsive molecules in each layer. We addedinterlayer regulatory connections between glucoseresponsive molecules in different layers using bioinformatics methods and information in public databases. The result was a glucose-responsive transomic network in the skeletal muscle of WT and ob/ob mice. We identified transomic subnetworks specific to WT mice (blue), ob/ob mice (red), and common to both mice (green). (B) Definition of glucose-responsive molecules using fold change and FDR-adjusted p value (q value). *q value < 0.1 and absolute log 2 fold change > 0.585. (C) Definition of T 1/2 , an index of the temporal rate of response, and rapid and slow glucose-responsive molecules using T 1/2 . This figure was modified from Figure 1   glucose-responsive metabolites from the skeletal muscles of WT and ob/ob mice following oral glucose administration (n = 5 mice per genotype at all time points). To investigate the changes from fasting state, two time courses for each metabolite were divided by the geometric mean of the values of fasted WT mice and ob/ob mice (0 min) and then log 2 -transformed. Metabolites were ordered by hierarchical clustering using Euclidean distance and Ward's method. Right: The bars in the heat map are colored according to the extent of glucose responsiveness, meaning the change from fasting state (0 min) in WT and ob/ob mice. Metabolites that showed an absolute log 2 fold change ≥ 0.585 (2 0.585 = 1.5) and an FDR-adjusted p value (q value) ≤ 0.1 at any time point (20,60, 120, 240 min) were defined as glucose-responsive: increased (orange), decreased (purple), or were unchanged (white). To define an increase or decrease in time courses with changes in both directions at different times, we used the direction of change compared to time 0 at the earliest time point that showed a significant change. Metabolites written in blue text indicate glucose-responsive metabolites specific to WT mice; red text, specific to ob/ob mice. (B) Increased and decreased metabolites in the skeletal muscles of WT mice and ob/ob mice. Blue, WT specific; red, ob/ob specific; gray, not glucose-responsive metabolites either in WT mice or in ob/ob mice. The number of each type of glucose-responsive metabolites and their percentages of the total quantified metabolites are shown. (C) Rapid, intermediate, and slow responses in increased metabolites specific to WT mice, decreased metabolites specific to WT mice, increased metabolites specific to ob/ob mice, and decreased metabolites specific to ob/ob mice. (D) Graphs showing the metabolites with responses that were specific to WT mice (blue boxes) and specific to ob/ob mice (red boxes). Blue lines are the responses of the WT mice and red lines are the responses of the ob/ob mice. Data are shown as the mean and standard error of the mean (SEM) of five mice per genotype. *q value < 0.1 and absolute log 2 fold change > 0.585. Metabolites are abbreviated as follows: F6P, fructose 6-phosphate; PEP, phosphoenolpyruvate; Glycerol-3P, glycerol 3-phosphate; 3-OH-butyrate, 3-hydroxybutanoate. www.nature.com/scientificreports/ Transcriptomics. To elucidate the transcriptional changes and controls in the skeletal muscle of WT and ob/ob mice after glucose administration, we performed transcriptomic analysis using RNA sequencing. Of the 14,978 genes analyzed, 4,264 that were significantly changed after oral glucose administration were identified as glucose-responsive genes (Fig. 3A,B; Data File S3). A heatmap of the glucose-responsive genes is shown in Fig. 3A. The responses were categorized into three groups (rapid, intermediate, or slow) according to their T 1/2 s as in the analysis of the metabolites (Fig. 3C,D). Pathway enrichment analysis was also performed for each type of response (Table 1 and Data File S4). We assigned glucose-responsive genes encoding metabolic enzymes to the Enzyme layer of the transomic network, and glucose-responsive genes encoding transcription factors to the TF layer of the transomic network ( Figs. 1 and 5). The number of upregulated and downregulated genes in WT and ob/ob mice is shown in Fig. 3B. The number of glucose-responsive genes specific to ob/ob mice (1414 upregulated, 995 downregulated) was larger than that specific to WT mice (624 upregulated, 472 downregulated). A total of 637 common genes were upregulated and 264 were downregulated in WT and ob/ob mice. The calculation of time constants revealed that the number of rapidly responding glucose-responsive genes was larger in ob/ob mice than in WT mice (Fig. 3C). Genes upregulated in both WT and ob/ob mice included those involved in central carbon metabolism, such as hexokinase 2 (Hk2), fatty acid synthase (Fasn), and stearoyl-coenzyme A (CoA) desaturase 1(Scd1), and the responses in ob/ob mice were larger than those in WT mice (Fig. 3E). Some genes involved in the insulin signaling pathway also showed upregulation common to both WT and ob/ob mice, such as insulin receptor substrate 1 (Irs1) and nitric oxide synthase 3 (Nos3) (Fig. 3E). Genes downregulated in both WT and ob/ob mice included those involved in oxidative phosphorylation such as NADH dehydrogenase (ubiquinone) iron-sulfur protein 4 (Ndufs4) (Fig. 3E). Genes specifically downregulated in WT mice contained pyruvate dehydrogenase kinase 4 (Pdk4) (Fig. 3E). Genes specifically upregulated in ob/ob mice were relatively enriched in pathways related to cell adhesion ( Table 1). The gene 3-hydroxybutyrate dehydrogenase 1 (Bdh1), which is involved in ketone body metabolism, was also specifically upregulated in ob/ob mice. Genes specifically downregulated in ob/ob mice included those involved in the TCA cycle such as succinyl-CoA synthetase beta subunit (Suclg2), and those involved in branched-chain amino acid (BCAA) degradation such as 2-oxoisovalerate dehydrogenase beta subunit (Bckdhb) (Fig. 3E). Genes specifically downregulated in ob/ob mice were relatively enriched in the proteasome pathway and ribosomal proteins (Table 1).
Next, we performed hierarchical clustering analysis of transcriptome data and bioinformatics analysis of the binding motifs of gene clusters using the transcription factor database TRANSFAC (Figs. 3E and S5A, B; Data Files S5 and S6) to estimate the regulatory connections between transcription factors and genes 15,16 . We predicted the regulatory connections between a transcription factor and a gene if the binding motifs of the transcription factor were enriched in the promoter regions of the genes in a cluster. For example, we inferred that early growth response protein 1 (Egr1) is a transcription factor that regulates some of the genes upregulated in WT and ob/ob mice (Fig. 3E). A comparison of the estimated regulatory connections with those predicted from chromatin immunoprecipitation (ChIP) experimental data from the ChIP-Atlas database (http:// chip-atlas. org/) 17 showed that the results from the two methods mostly overlapped ( Fig. S5C; Data File S7). The estimated regulatory connections between the transcription factors and the genes encoding metabolic enzymes acted as connections between the TF layer and the Enzyme layer in the transomic network.

Phosphorylation of insulin signaling molecules.
Phosphorylation is an important factor for regulating metabolic reactions. Direct phosphorylation of an enzyme can regulate its activity, and phosphorylation of a transcription factor can regulate the expression level of downstream enzymes. Therefore, we measured the phosphorylation of 10 enzymes, transcription factors, and signaling molecules in the insulin pathway by performing western blot analysis of protein samples prepared from the skeletal muscle of WT and ob/ob mice during oral glucose administration ( Fig. S6; Data File S8). The band intensities were quantified, and the results were used to determine if the phosphorylation was glucose-responsive.
We were able to detect many glucose-responsive phosphorylated proteins from the analysis (Fig. 4). The level of phosphorylated ribosomal protein S6 was increased in both WT and ob/ob mice. Phosphorylated Akt was specifically increased in WT mice, and phosphorylated glycogen phosphorylase was specifically decreased in WT mice. Glycogen synthase kinase 3 β (Gsk3β) and cAMP response element-binding protein (Creb) were specifically increased in ob/ob mice. Some molecules showed the opposite responses in WT and ob/ob mice. For example, phosphorylated forkhead box protein 1 (Foxo1) was transiently increased in WT mice but decreased in ob/ob mice; phosphorylated glycogen synthase (Gs) was decreased in WT mice and increased in ob/ob mice. Phosphorylated extracellular signal-related kinase (Erk) and AMP-activated protein kinase α (Ampkα) was not affected by glucose administration in both WT and ob/ob mice. In the subsequent transomic analysis, metabolic enzymes with glucose-responsive phosphorylation were assigned to the Enzyme layer, and transcription factors with glucose-responsive phosphorylation were assigned to the TF layer.
Regulatory glucose-responsive transomic network. A regulatory transomic network of glucoseresponsive molecules in the skeletal muscle was constructed with five layers: Insulin signal, TF, Enzyme, Reaction, and Metabolite ( Fig. 5; Data File S9). We constructed the transomic network in the skeletal muscle using a method we previously developed for the transomic network in the liver 12 . Briefly, glucose-responsive molecules were assigned to the corresponding layers as nodes, and the edges between the nodes were drawn to show the interlayer regulatory connections of glucose-responsive molecules retrieved from pathway databases such as Kyoto Encyclopedia of Genes and Genomes (KEGG) and Braunschweig Enzyme Database (BRENDA) [18][19][20] (Fig. 5A). www.nature.com/scientificreports/ By constructing regulatory transomic networks in WT and ob/ob mice, we were able to identify WT specific, ob/ob specific, and common responses of molecules and interlayer regulatory connections to glucose administration ( Fig. 5B; green, common; blue, WT specific; red, ob/ob specific). In the Metabolite layer, the number of WT mice specific glucose-responsive molecules was larger than ob/ob mice specific glucose-responsive molecules, and no molecules responded commonly in WT and ob/ob mice. Therefore, most of the interlayer regulatory connections between the Metabolite layer and the Reaction layer were specific to WT mice, suggesting that metabolic regulation by a metabolite itself after glucose administration is impaired in obesity. By contrast, approximately 55% of glucose-responsive genes in the Enzyme layer and the interlayer regulatory connections between the Enzyme layer and the Reaction layer were classified as ob/ob specific, suggesting that transcriptional regulation compensated for the regulation by metabolites that was lost in obese mice. The number of common glucose-responsive genes in the Enzyme layer and its regulatory connections was approximately 40% of the ob/ob specific ones.
The numbers of glucose-responsive metabolic reactions regulated by metabolites (Metabolite layer), genes (Enzyme layer), or both were calculated (Fig. 5C). The results suggested that the metabolic reactions in WT mice were mainly regulated by metabolites, and those in ob/ob mice were mainly regulated through gene expression. We also classified the regulators of metabolic reactions according to their time constants (T 1/2 ), and revealed that a large number of metabolic reactions was affected by the rapidly responding (< 20 min) metabolites and genes in both the WT and ob/ob networks (Fig. 5D). Glucose-responsive metabolites specific to WT mice included cofactors such as ATP, AMP, and UDP, which could have a large effect on the Reaction layer (Fig. 5E).

Comparison of the regulatory transomic networks of WT and ob/ob mice.
To analyze how each metabolic pathway was regulated in the regulatory transomic networks of WT and ob/ob mice, we constructed a simplified transomic network using a method that we previously developed 12 ( Fig. 6; Data File S10). Briefly, we converted the Reaction layer into the Pathway layer by placing metabolic reactions in a specific metabolic pathway into a single metabolic pathway node, according to the KEGG metabolic pathway.
In WT mice, various metabolic pathways were regulated by metabolites (Fig. 6A). In particular, carbohydrate metabolic pathways were regulated by WT specific glucose-responsive metabolites such as ATP, citrate, and phosphoenolpyruvate (PEP). Although the effects of glucose-responsive genes encoding metabolic enzymes were smaller than the metabolites, some lipid metabolic pathways such as glycerolipid and glycerophospholipid metabolisms were more strongly regulated by glucose-responsive genes than others (Fig. S7A). In ob/ob mice, the regulation of glucose-responsive metabolites was decreased and that of glucose-responsive genes encoding metabolic enzymes was increased (Fig. 6B). The decreased regulation by metabolites was particularly large in carbohydrate metabolic pathways (Fig. S7A). Regulation by glucose-responsive genes was increased in most carbohydrate and lipid metabolic pathways, with the exception of glycerolipid metabolism. Amino acid metabolic pathways showed relatively small changes in the percentage of metabolic reactions regulated by glucoseresponsive metabolites and genes.
Glycolysis, TCA cycle, BCAA degradation, and ketone body metabolism. Finally, we focused on metabolic pathways and their regulatory networks related to glucose (Fig. 7). abundance for 4,264 glucose-responsive genes in the skeletal muscles of WT and ob/ob mice following oral glucose administration (n = 11 or 12 mice per genotype at 0 min, n = 5 mice per genotype at 20 min, n = 5 mice per genotype at 60 min, n = 5 mice per genotype at 120 min, and n = 5 mice per genotype at 240 min). To investigate the changes from fasting state, two time courses for each gene were divided by the geometric mean of the values of WT mice and ob/ob mice in the fasting state (0 min), and then log 2 -transformed. Genes were ordered by hierarchical clustering using Euclidean distance and Ward's method (Data File S3). Right: The bars in the heat map are colored according to glucose responsiveness, meaning the change from fasting state (0 min) in WT and ob/ob mice. Genes that showed an absolute log 2 fold change ≥ 0.585 (2 0.585 = 1.5) and a q value ≤ 0.1 at any time point (20, 60, 120, 240 min) were defined as glucose-responsive: upregulated (orange), downregulated (purple), or were unchanged (white). (B) Increased and decreased genes in the skeletal muscle of WT mice and ob/ob mice. Blue, WT specific; red, ob/ob specific; green, glucose-responsive genes common to both; gray, not glucose-responsive metabolites either in WT mice or in ob/ob mice. The number of each type of glucose-responsive genes and their percentages out of the total quantified genes are shown. (C) Rapid, intermediate, and slow responses in increased genes specific to WT mice, decreased genes specific to WT mice, increased genes specific to ob/ob mice, and decreased genes specific to ob/ob mice. (D) Rapid, intermediate, and slow responses in increased and decreased genes common to both WT and ob/ob mice. (E) Graphs showing the gene expression time courses for the indicated genes. Genes include those that exhibited changes in common to both WT and ob/ob (green boxes), changes specific to WT mice (blue box), and changes specific to ob/ob mice (red boxes). Within the graphs, blue lines are the responses of the WT mice and red lines are the responses of the ob/ob mice. The inferred regulatory connections are shown as arrows from transcription factors to genes. The regulatory connections were inferred using hierarchical clustering analysis of gene expression time courses together with a transcription factor database TRANSFAC 15,16 . See Figure S5 for statistical confidence in the inferred transcription factors; see Data File S5 for the unabbreviated names of the transcription factors. Data are shown as the mean and SEM (n = 11 or 12 mice per genotype at 0 min, n = 5 mice per genotype at 20 min, n = 5 mice per genotype at 60 min, n = 5 mice per genotype at 120 min, and n = 5 mice per genotype at 240 min). *q value < 0.1 and absolute log 2 fold change > 0.585. www.nature.com/scientificreports/ Glycolysis. In WT mice, although blood glucose levels increased after glucose administration, most metabolites in glycolysis were not defined as "glucose-responsive. " The glycolysis network contained many allosteric inhibitors that decreased specifically in WT mice, such as ATP and citrate (Fig. 7A). We also found upregulation in some glycolytic genes such as Hk2, and downregulation in Pdk4, which inhibits pyruvate dehydrogenase by phosphorylation 21 .
In ob/ob mice, F6P and lactate were defined as specifically increased metabolites, and some gene expression such as Hk2 showed a larger increase than in WT mice. Most of the responses of allosteric regulators of glycolysis were lost in ob/ob mice. . Glucose 6-phosphate (G6P) was not defined as a glucose-responsive molecule (q value at 60 min = 0.14), but its time series was highly correlated with F6P (Pearson's r = 0.99).
TCA cycle. In WT mice, four metabolites in the TCA cycle, namely citrate, cis-aconitate, succinate, and malate, decreased after oral glucose administration (Fig. 7B). Although fumarate was not defined as a glucose-responsive molecule (q value at 60 min = 0.13), its time series was highly correlated with malate (Pearson's r = 0.96) [21][22][23][24] . In ob/ob mice, the abundance of some metabolites was smaller than that in WT mice before glucose administration, and the metabolites did not show a large response to glucose. Some studies have reported a decrease in intermediates of the TCA cycle in the skeletal muscle of obese mice 25,26 .
BCAA degradation. BCAA degradation pathway and its regulatory network included some glucose-responsive molecules in ob/ob mice (Fig. 7C). Valine showed a rapid decrease after oral glucose administration. Leucine and isoleucine were not defined as glucose-responsive molecules (q value at 20 min = 0.15, 0.16), but their time series were highly correlated with valine (Pearson's r = 0.98 for leucine, 0.98 for isoleucine). Some genes involved in BCAA degradation, such as Bckdhb, showed a rapid downregulation. Bckdh kinase (Bckdk) inhibits Bckdh by phosphorylation 27 , which was not defined as a glucose-responsive molecule (q value at 60 min = 0.11), but its time series was negatively correlated with Bckdhb expression (Pearson's r = -0.96). We found a similar decrease in Suclg2 in the TCA cycle, which metabolizes succinyl CoA, one of the BCAA degradation products (Fig. 7B). In WT mice, BCAAs were not defined as glucose-responsive molecules (q value at 20 min = 0.15 to 0.23), but their time series showed a positive correlation with those in ob/ob mice (Pearson's r = 0.77 to 0.90). Acyl-CoA dehydrogenase short/branched chain (Acadsb), methylcrotonoyl-CoA carboxylase 2 (Mccc2), and Suclg2 showed significant decrease in WT mice, but only Acadsb was defined as glucose-responsive because of its fold changes (fold change at 20 min = 0.64, 0.76, 0.80).
Ketone body metabolism. In WT mice, 3-OH butyrate, a ketone body, showed a rapid and strong decrease (0.13-fold at 20 min) (Fig. 7D). The network included the decreases in metabolites in the TCA cycle as allosteric inhibitors of the metabolic enzyme that degrades acetoacetate. In ob/ob mice, 3-OH butyrate did not show a significant decrease (q value at 60 min = 0.12), but Bdh1 was rapidly upregulated. www.nature.com/scientificreports/

Discussion
In this study, we performed transomic analysis of the skeletal muscles obtained from WT and ob/ob mice after the oral glucose tolerance test to construct a large-scale glucose-responsive regulatory network of metabolism. In WT mice, the number of glucose-responsive metabolites was about 2.5-fold larger than that in ob/ob mice, and many metabolic reactions were affected by these glucose-responsive metabolites. In particular, the responses of cofactors such as ATP, and TCA cycle intermediates such as citrate and succinate, might affect carbohydrate and amino acid metabolism. By contrast, the number of glucose-responsive genes encoding metabolic enzymes in ob/ob mice was about 1.8-fold larger than that in WT mice, and the genes were mainly related to carbohydrate and lipid metabolism. We also found some characteristic glucose-responsive regulatory pathways in central carbon, branched amino acids, and ketone body metabolism. The WT mice showed few significant changes in the metabolites of glycolysis despite the administration of glucose. This can be attributed to two hypotheses. The first hypothesis is that the influx into the muscle glycolytic system does not increase as much as the blood glucose levels. A recent study showed that the influx of orally administered labelled glucose into the glycolysis of gastrocnemius muscle (white muscle), which was used in this study, is much smaller than that of soleus muscle (red muscle) 28 . The second hypothesis is that the efflux of glycolysis increased by activation as much as the influx. Thus, the network of glycolysis in WT mice included some activating regulatory pathways, such as the decrease of allosteric inhibitor including ATP and TCA cycle intermediates, the upregulation of glycolytic genes, and the downregulation of Pdk4 (Fig. 7A), which could contribute to blood glucose homeostasis after glucose administration. Because blood lactate increased (Fig. S4B), much of the increased glucose flowed into glycolysis after administration might be released into the blood as lactate 29,30 . This is consistent with the decrease in ATP and TCA cycle intermediates, which suggests a decrease in TCA cycle flux (Fig. 7B). By contrast, ob/ob mice showed the considerable increase in F6P and G6P, suggesting that the conversion from glucose to G6P was activated in response to blood glucose levels, but to a lesser extent in the glycolytic system. Consistently, in ob/ob mice, most of allosteric activation of glycolysis was lost, and Hk2 showed larger increase than WT mice, which could be compensational regulation www.nature.com/scientificreports/ for allosteric regulation. In this study, some amino acids including BCAA in the blood and skeletal muscle were decreased after glucose administration similar to the effect on the liver 12 , suggesting suppression of protein degradation and promotion of protein synthesis in the insulin target organs 22,22,24,31 . In BCAA degradation pathway, as well as the decrease of the amino acids, we found the transcriptional repression of the metabolic enzymes, including Bckdhb, and the transcriptional activation of Bckdk, a inhibitory kinase of Bckdh 27 , in ob/ob mice and partially in WT mice (Fig. 7C). These responses might suppress the degradation of BCAAs to the TCA cycle intermediates and contribute to the decrease in TCA cycle flux in the skeletal muscle. The blood level of a ketone body, an alternative energy source in the fasting state, was decreased in both WT and ob/ob mice after glucose administration. We also found that ketone levels in skeletal muscle showed a similar time series as those in the blood (Fig. 7D), suggesting that intramuscular degradation of ketone body to acetyl CoA was also reduced. Decreased degradation of these metabolites could contribute to a decrease in TCA cycle intermediates, but The regulatory transomic network for glucose-responsive metabolic reactions. The left diagram represents the network as colored nodes in the layers and edges between the layers with colored nodes representing glucose-responsive molecules and colored edges representing interlayer regulatory connections: green, glucose-responsive molecules and interlayer regulation common in both WT and ob/ob mice; blue, specific to WT mice; red, specific to ob/ob mice; pink, opposite responses between WT and ob/ob mice. The numbers of each type of glucose-responsive node and edge are shown in the same colors in the network summary to the right. The Insulin signal layer was the insulin signaling pathway constructed in our previous phosphoproteomic study 11 . The Enzyme, Reaction, and Metabolite layers were organized into global metabolic pathways (mmu01100) in the KEGG database 18,19 . www.nature.com/scientificreports/ further research is needed to understand why the reduction was specific to WT mice. As mentioned above, the WT-specific responses of ATP, a product of the TCA cycle, and TCA cycle intermediates could activate glycolysis and glucose metabolism. Thus, elucidating the mechanism of the response will be helpful in understanding blood glucose homeostasis and its disruption by obesity. We previously constructed a glucose-responsive transomic network in the liver of WT and ob/ob mice 12 . The liver network contained more glucose-responsive molecules and regulatory connections than the skeletal muscle network, but the differences between WT and ob/ob mice were similar between the liver and skeletal muscle. In both organs, many metabolic reactions in the WT networks were regulated by metabolites, whereas in the ob/ob networks, much of the regulation by metabolites was lost and metabolic regulation by gene expression was activated. There were also similarities in the regulation of the metabolic pathway, such as the regulation of carbohydrate metabolism by metabolites and the regulation of lipid metabolism by gene expression. We are currently performing a detailed comparative analysis between the liver network and skeletal muscle network.
To construct a comprehensive glucose-responsive network, it was necessary to integrate more omics data into our network. Because the Insulin signal layer was determined by western blot analysis, the numbers of glucoseresponsive molecules and regulatory connections of the layer were very limited compared to those of the other layers. Integration of phosphoproteomic data and kinase-substrate interactions will facilitate a more extensive The condensed regulatory transomic network of the response to glucose in WT and ob/ob mice. The color of nodes (glucose-responsive molecules) and edges (interlayer regulatory connections) indicate the type of molecules and regulation as described in Fig. 5B. In the TF layer, the transcription factors encoded by glucose-responsive genes are italicized, and those showing glucose-responsive phosphorylation have the prefix "p. " The Enzyme layer contained only the metabolic enzymes that were regulated by glucose-responsive changes in gene expression, and not those regulated only by phosphorylation. Two types of metabolic pathway nodes were included: the pathway that exhibited significant associations with any glucose-responsive molecule (Fig. S7B); and the pathway whose percentage of regulated reactions was in the top 10% by either glucose-responsive metabolites or by glucose-responsive genes encoding metabolic enzymes. Dashed boxes enclose the nodes for the lipid, carbohydrate, and amino acid classes. Glucose-responsive metabolites that exhibited significant associations with any metabolic pathway were included. Glucose-responsive transcription factors that regulated five or more metabolic reactions were included. The interlayer regulatory connections from the Metabolite to the Pathway layer included only those that regulate five or more metabolic reactions. The size of the nodes and width of the edges indicate the relative number of the regulated metabolic reactions. See Data File S10 for the unabbreviated names of the metabolic pathway nodes. www.nature.com/scientificreports/ evaluation of the effects from the Insulin signal layer to the Reaction layer [32][33][34] . The transcription factors of the glucose-responsive genes were determined based on the binding motifs in the promoter sequences and the temporal patterns. Because not all motifs are bound by transcription factors, direct measurements of transcription factor binding using ChIP sequencing analysis will identify a more accurate and extensive regulatory network of glucose-responsive genes 17,35,36 . Furthermore, to accurately identify the effects of obesity on the glucoseresponsive network, it is necessary to compare the glucose-responsive networks from various mouse models of obesity, such as HFD mice, as well as ob/ob mice. Although our transomic network was not comprehensive, we revealed several important features of metabolic regulation in the skeletal muscle after glucose administration. An extension of this in vivo transomic analysis will lead to a better understanding of glucose homeostasis at the whole-body level and its dysregulation in obesity.

Materials and methods
Animals and sample preparation. Animal experiments were performed as previously described 12 .
C57BL/6J WT mice or ob/ob mice at ten weeks of age were purchased from Japan SLC Inc. (Shizuoka, Japan). The phenotypic data of the mice are summarized in Table S1. Animal experiments were approved by the animal ethics committee of The University of Tokyo and according to the ARRIVE guidelines and the University of Tokyo guidelines for the care and use of laboratory animals. After purchasing, the mice were housed 2-3 to a cage and fasted overnight at 23 °C (6 p.m. to 10 a.m.). Overnight-fasted mice were administered an oral glucose load of 2 g/kg body weight (41.8-47.8 mg for WT mice, 76.8-94.4 mg for ob/ob mice). To measure blood glucose and insulin levels, 15 μL blood was collected from the tail veins at 0, 2, 5, 10, 15, 20, 30, 45, 60, 90, 120, 180, and 240 min after glucose administration (n = 5). We used the blood glucose and insulin levels measured in our previous study 12 (Fig. S2). For the metabolome and transcriptome studies, mice were sacrificed at 0, 20, 60, 120, and 240 min after glucose administration, and the gastrocnemius muscle was excised. Muscle samples were frozen immediately in liquid nitrogen and homogenized with dry ice. The powdered samples were divided and used for metabolomics, lipidomics, transcriptomics, a glycogen assay, and western blotting.
Metabolomics. Metabolomic analysis was performed as previously described 12   . Regulatory transomic network for glucose-responsive metabolic reactions in glycolysis, TCA cycle, BCAA degradation, and ketone body metabolism. The regulatory transomic network for glucose-responsive metabolic reactions in glycolysis (A), TCA cycle (B), BCAA degradation (C), and ketone body metabolism (D) in the skeletal muscle of WT mice and ob/ob mice. The information for the pathways was obtained from "glycolysis/gluconeogenesis" (mmu00010), "citrate cycle (TCA cycle)" (mmu00020), "valine, leucine and isoleucine degradation" (mmu00280), and "synthesis and degradation of ketone bodies" (mmu00072) in the KEGG database 18,19 . Graphs of the time courses of measured molecules are shown for corresponding nodes as the means and SEMs (n = 5 mice per genotype for metabolite, n = 11 or 12 mice per genotype for gene expression at 0 min, n = 5 mice per genotype for gene expression at 20, 60, 120, 240 min, n = 5 mice per genotype for phosphorylation). The colors of the frames indicate WT mice-specific glucose-responsive molecules (blue), ob/ob mice-specific glucose-responsive molecules (red), and common glucose-responsive molecules (green). The dashed frames indicate molecules that were not included in the glucose-responsive transomic network. Diamond nodes indicate metabolic reactions. The colored edges indicate interlayer regulatory connections: WT mice-specific regulatory connections (blue), ob/ob mice-specific regulatory connections (red), and common interlayer regulatory connections (green). From the Metabolite to Reaction layers, only allosteric regulatory connections are colored. Black edges indicate the relationship between metabolic reactions and its substrate/ product. The reversibility of metabolic reactions was obtained from the KEGG database 18,19  www.nature.com/scientificreports/ by normalization to the reference compound 3-aminopyrrolidine. The metabolites were identified by comparing their m/z values and relative migration times to the metabolite standards. Quantification was performed by comparing peak areas to calibration curves generated using internal standardization techniques with methionine sulfone. The other conditions were identical to those previously described 38 . To analyze anionic metabolites, a commercially available COSMO(+) (chemically coated with cationic polymer) capillary (50 µm i.d. × 105 cm) (Nacalai Tesque, Kyoto, Japan) was used with a 50 mM ammonium acetate solution (pH 8.5) as the electrolyte. Methanol/5 mM ammonium acetate (50% v/v) containing 0.01 µM hexakis(2,2-difluoroethoxy)phosphazene was delivered as the sheath liquid at 10 µL/min. ESI-TOF MS was performed in the negative ion mode, and the capillary voltage was set to 3.5 kV. For anion analysis, trimesate and D-camphor-10-sulfonic acid were used as the reference and internal standard, respectively. The other conditions were identical to those described previously 39 . Agilent MassHunter software (Agilent technologies) was used for data analysis [38][39][40] . We used the blood metabolome data obtained in our previous study 12 . Lipidomics. Lipidomic analysis was performed as previously described 13 . Lipidomic profiling of the skeletal muscle was performed by Metabolon, Inc. Lipids were extracted from samples with dichloromethane and methanol using the modified Bligh and Dyer procedure in the presence of internal standards, with the lower organic phase used for analysis. The extracts were concentrated under nitrogen and reconstituted in 0.25 mL dichloromethane:methanol (50:50) containing 10 mM ammonium acetate. The extracts were placed in vials for infusion-MS analyses, which were performed on the SelexION equipped Sciex 5500 QTRAP mass spectrometer using both the positive and negative ion modes. Each sample was subjected to two analyses, with ion mobility spectrometry-MS conditions optimized for lipid classes monitored in each analysis. The 5500 QTRAP was operated in the multiple reaction monitoring mode to monitor the transitions for more than 1,100 lipids from up to 14 lipid classes. Individual lipid species were quantified based on the ratio of the signal intensity for target compounds to the signal intensity for an assigned internal standard of known concentration. Fourteen lipid class concentrations were calculated from the sum of all molecular species within a class.
Glycogen assay. Glycogen content was determined as previously described with some modifications 41 .
Approximately 20 mg of the skeletal muscle was digested with 1.2 mL of 30% (w/v) potassium hydroxide solution for 1 h at 95 °C and neutralized with 61.2 μL glacial acetic acid. The total protein concentration of the muscle digest was determined by the BCA assay and adjusted to 1 μg protein/μL. Glycogen was extracted from the digested skeletal muscle using Bligh and Dyer method to remove lipids 42 . The digested skeletal muscle (50 μL) was mixed with 120 μL ice-cold methanol, 50 μL chloroform, 10 μL of 1% (w/v) linear polyacrylamide, and 70 μL water. After incubation on ice for 30 min, the mixture was centrifuged at 12,000×g to remove the separated aqueous layer. The glycogen was precipitated by the addition of 200 μL methanol and centrifugation at 12,000×g for 30 min at 4 °C, washed with ice-cold 80% (v/v) methanol, and dried completely. Glycogen pellets were suspended in 20 μL of 0.1 mg/mL amyloglucosidase (Sigma-Aldrich) in 50 mM sodium acetate buffer and incubated for 2 h at 55 °C to digest the glycogen. The concentration of the glucose produced from the glycogen was determined using the Amplex Red Glucose/Glucose Oxidase Assay kit (Thermo Fisher Scientific), according to the manufacturer's instructions.
Transcriptomics. Transcriptomic analysis was performed as previously described 12 . Total RNA was extracted from the skeletal muscle using the RNeasy Mini Kit (QIAGEN) and QIAshredder (QIAGEN); the quantity was assessed using the Nanodrop (Thermo Fisher Scientific) and the quality was assessed using the 2100 Bioanalyzer (Agilent Technologies). cDNA libraries were prepared using the SureSelect strand-specific RNA library preparation kit (Agilent Technologies). The resulting cDNAs were subjected to 100 base paired-end sequencing on the Illumina HiSeq2500 Platform (Illumina) 43 . Sequences were aligned to the mouse reference genome obtained from the Ensembl database 44,45 (GRCm38/mm10, Ensembl release 97) using the STAR software package TopHat (v.2.5.3a) 46 . The RSEM tool (v.1.3.0) was used to assemble transcript models (Ensembl release 97) from aligned sequences and to estimate gene expression level 47 . Gene expression level was shown as fragments per kilobase of exon per million mapped fragments.
Western blot analysis. Total proteins were extracted from the skeletal muscle with methanol:chloroform:water (2.5:2.5:1). Ice-cold methanol was added to the skeletal muscle at a concentration of 100 mg/mL of the weight of the skeletal muscle, and the suspension (400 μL) was mixed with chloroform (400 μL) and water (160 μL), followed by centrifugation at 4600×g for 10 min at 4 °C. The aqueous and organic phases were removed and 800 μL ice-cold methanol was added to the interphase to precipitate proteins. The resulting pellet was suspended with 400 μL lysis buffer ( www.nature.com/scientificreports/ Identification of glucose-responsive molecules. Glucose-responsive molecules were determined as previously described 12 . Molecules that were detected in less than half of the replicates in either WT or ob/ob mice at any time point after oral glucose administration were removed from the analysis. A molecule with a statistically significant change in response to oral glucose administration was defined as a glucose-responsive molecule according to the following criteria. The fold change of the mean amount at each time point over the mean amount at fasting state (0 min) was calculated for each molecule. The significance of change at each time point was tested by the two-tailed Welch's t-test for each metabolite and phosphorylation, and by the edgeR package (version 3.26.8) of the R language (version 3.6.1) with the default parameters for each gene 48 . Metabolite, gene, and phosphorylation that showed an absolute log 2 fold change ≥ 0.585 (2 0.585 = 1.5) and an FDR-adjusted p value (q value) ≤ 0.1 at any time point were defined as a glucose-responsive metabolite ( Fig. 2A,B), gene (Fig. 3A,B), and phosphorylation (Fig. 4). The q values were calculated by Storey's procedure 49 . To define an increase or decrease in time courses with changes in both directions at different times, we used the direction of change compared to time 0 at the earliest time point that showed a significant change.
Clustering analysis. Time courses for each metabolite of WT mice and ob/ob mice were normalized by dividing by the geometric mean of the values of WT mice and ob/ob mice in the fasting state (0 min) followed by log 2 transformation. We combined the two time courses of WT and ob/ob mice for each metabolite and performed hierarchical clustering of the combined time courses using Euclidean distance and Ward's method (Fig. S3). Based on the clustering tree, we defined eight different clusters of metabolites, showing similar or different responses between WT and ob/ob mice. Clustering analysis of gene expression was performed as previously described with some modifications 12 . Time courses for the expression of each gene of WT and ob/ob mice were normalized by subtracting the average expression values of the time courses of both mice and then dividing the resulting values by the standard deviation (Z-score normalization). We combined the two time courses of WT and ob/ob mice for each gene and performed hierarchical clustering of the combined time courses using Euclidean distance and Ward's method (Fig. S7A). The genes with significant differences between WT and ob/ob mice before glucose administration (0 min) (q value < 0.1) or a significant response at any time point in either WT or ob/ob mice (q value < 0.1) were selected for the clustering analysis (12,301 genes). For the selection, the p value was calculated using the edgeR package (version 3.26.8) of the R language (version 3.6.1) with the default parameters 48 , and the q value was calculated by Storey's procedure 49 . Pathway enrichment analysis. We performed pathway enrichment analysis of glucose-responsive genes (Table 1; Data File S4). The enrichment of the genes in each pathway was determined using the one-tailed Fisher's exact test. We used the genes detected in more than half of the replicates in WT and ob/ob mice at all time points as background. We used the pathways in Metabolism, Genetic Information Processing, and Cellular Processes from the KEGG database 18,19 . Prediction of the transcription factor binding motif and inference of regulatory connections between transcription factors and genes. Analysis of transcription factors was performed as previously described 12 . The flanking regions around the major transcription start site of genes were extracted from GRCm38/mm10 (Ensembl, release 97) using Ensembl BioMart 50 . The region from -300 bp to + 100 bp of the major transcription start site was defined as the flanking region, according to FANTOM5 analysis of the time course 51 . The transcription factor binding motifs in each flanking region (Fig. S5B) were predicted using TRANS-FAC Pro, a transcription factor database, and Match, a transcription factor binding motif prediction tool 15,16 . The threshold for each transcription factor binding motif prediction was set using extended vertebrate_non_redun-dant_min_FP.prf, a parameter set in TRANSFAC Pro 12 .
For the inference of regulatory connections between transcription factors and genes, we performed transcription factor motif enrichment analysis of the genes in each cluster (Fig. S5B). The enrichment of transcription factor binding motif in the flanking regions of genes in each cluster was determined by the one-tailed Fisher's exact test, and transcription factor binding motifs with q value ≤ 0.1 were defined as significantly enriched. The q values were calculated by the Benjamini-Hochberg procedure 52 . We used the genes analyzed in the hierarchical clustering as background. To reduce the number of statistical tests, the clusters that contained ≥ 100 genes were analyzed. If a transcription factor binding motif was enriched in the promoter regions of the genes in a cluster, we inferred the regulatory connections between the corresponding transcription factor and the genes in the cluster. To avoid overestimation, we excluded a cluster from the inference if the transcription factor binding motif was more enriched in the children clusters that contained ≥ 100 genes. To compare the enrichment of transcription factor binding motifs between clusters, we calculated the odds ratio of the transcription factor binding motifs for each cluster.
For validation of the inferred regulatory connections, we examined the overlap between the inferred genes of each transcription factor and those predicted from experimental ChIP data from the ChIP-Atlas database 17 (Fig. S5C). The genes for which ChIP sequencing peaks of a transcription factor were detected in the flanking region around the transcription start sites were obtained using "Target Genes, " a prediction tool in the ChIP-Atlas. We used the flanking regions from − 1000 to + 1000 bp of the transcription start sites in Target Genes. The overlap between the inferred genes and genes from ChIP data was determined by the one-tailed Fisher's exact test, and those with q value ≤ 0.1 were defined as significant. The q values were calculated by the Benjamini-Hochberg procedure 52  www.nature.com/scientificreports/ Insulin signaling pathway. The insulin signaling pathway in Fig. 4 is a subset of the nodes of the insulin signaling pathway in the KEGG database (mmu04910) 18,19 . We added regulatory input to Creb from the PI3K-Akt signaling pathway (mmu04151), MAPK signaling pathway (mmu04010), and AMPK signaling pathway (mmu04152), and regulatory input to FoxO1 from the FoxO signaling pathway (mmu04068) in the KEGG database. The edges from Akt to Ampk and from p38 to insulin receptor substrate were added according to previous studies 53,54 .
Construction of the regulatory glucose-responsive transomic network. The transomic network was constructed as previously described with some modifications 12 . The regulatory glucose-responsive transomic networks consisted of five layers, namely Insulin signal, TF, Enzyme, Reaction, and Metabolite, with interlayer regulatory connections (Fig. 5A,B). The Insulin signal layer is the insulin signaling pathway constructed in our previous phosphoproteomic study 11 . We included in the Insulin signal layer signaling molecules that we analyzed by western blotting; we did not include transcription factors such as Foxo1, or metabolic enzymes such as Gs in this layer. The TF layer consisted of all transcription factors with an inferred regulatory connection (Fig. S7B). The Enzyme layer consisted of all metabolic enzymes in the pathways in Metabolism obtained from the KEGG database 18,19 . The Reaction layer consisted of the metabolic reactions (based on EC number) corresponding to the metabolic enzymes in the Enzyme layer. The Metabolite layer consisted of all metabolites analyzed by CE-MS. Only the molecules and reactions corresponding to genes that were expressed in at least one sample were included in the Insulin signal, TF, Enzyme, and Reaction layers. Not all 15,608 genes were included in the network. Glucose-responsive molecules were assigned to the corresponding layers as nodes. The Insulin signal layer consisted of insulin signaling molecules with glucose-responsive phosphorylation. The TF layer consisted of transcription factors encoded by glucose-responsive genes or those with glucose-responsive phosphorylation. The Enzyme layer consisted of metabolic enzymes encoded by glucose-responsive genes or those with glucoseresponsive phosphorylation. The Reaction layer consisted of "glucose-responsive metabolic reactions, " which were defined as metabolic reactions regulated by glucose-responsive molecules. The Metabolite layer consisted of glucose-responsive metabolites. We also determined the direction of glucose responsiveness. To determine a direction for time courses with both increased and decreased time points, we used the direction of change at the earliest time point with a significant difference from time 0 (fasting state). We did not determine a direction (increase or decrease) for metabolic reactions because we did not measure metabolic reaction activity.
To determine regulatory connections from the Enzyme and Metabolite layers to the Reaction layer, both the target of the regulatory connection (a metabolic reaction) and the regulating molecule (enzyme or metabolite) had to be glucose-responsive. Among the Insulin signal, TF, and Enzyme layers, the interlayer regulatory connections were determined using the directions of glucose responsiveness of the regulating molecule and the regulated molecules, and the types of interlayer regulatory connections, which were designated as either positive or negative. We defined positive interlayer regulatory connections as when both the regulating molecule and regulated molecule showed the same direction of change, namely, both increased or both decreased. We defined negative interlayer regulatory connections as when the regulating molecule and regulated molecule showed responses in the opposite direction, namely, one increased and the other decreased.
The interlayer regulatory connections between glucose-responsive molecules were determined according to databases. The interlayer connections from the Insulin signal layer to the TF layer were determined by the regulation of transcription factors by kinases retrieved from the KEGG database 18,19 . The interlayer connections from the TF layer to the Enzyme layer were determined from inferred regulatory connections between transcription factors and genes (Fig. 3E). The interlayer connections from the Enzyme layer to the Reaction layer were determined by connecting metabolic reactions to their corresponding metabolic enzymes according to the KEGG database 18,19 . The interlayer connections from the Metabolite layer to the Reaction layer comprised two types of regulatory connections: those mediated by allosteric regulators, which were retrieved from the BRENDA database 20 , and those mediated by the substrate or product of the reaction, which were retrieved from the KEGG database 18,19 . The types of regulatory connections made by glucose-responsive transcription factors were defined according to the Gene Ontology (GO) annotations obtained from the Mouse Genome Database 55 (Data File S6). The transcription factors that were included in the list of DNA-binding transcription repressors (GO:0001227) and not in the list of DNA-binding transcription activators (GO:0001228) were defined as transcription repressors. Foxo1 was added to the list of transcription activators based on previous studies of gluconeogenesis 56,57 . The effects of the phosphorylation of transcription factors on the types of regulatory connections were defined according to the KEGG database 18,19 . We used the allosteric regulation reported for mammals (Bos taurus, Felis catus, Homo sapiens, "Macaca, " "Mammalia, " "Monkey, " Mus booduga, Mus musculus, Rattus norvegicus, Rattus rattus, Rattus sp., Sus scrofa, "dolphin, " and "hamster") according to the BRENDA database 20 . Because the reversibility of metabolic reactions was not determined, metabolic reactions were assumed to be regulated by both the substrate and product.
Generation of a condensed transomic network based on metabolic pathway information. We condensed the regulatory transomic networks as previously described with some modifications 12 . First, we grouped the related metabolic reactions in a specific metabolic pathway into one "metabolic pathway node" (Pathway layer), and classified the metabolic pathway nodes into three classes-carbohydrate, lipid, and amino acid-according to the KEGG database 18,19 . Second, we selected two types of metabolic pathway nodes: one was a pathway that exhibited significant associations with any glucose-responsive metabolites or transcription factors; the other was a pathway whose percentage of regulated reactions was in the top 10% either by glucoseresponsive metabolites or by glucose-responsive genes encoding metabolic enzymes (Fig. S7A) www.nature.com/scientificreports/ between the metabolic reactions in a metabolic pathway and those regulated by a glucose-responsive molecule was determined by the one-tailed Fisher's exact test, and associations with a q value ≤ 0.01 were defined as significant. The q values were calculated by the Benjamini-Hochberg procedure 52 . We also selected glucose-responsive metabolites that exhibited significant associations with any metabolic pathway nodes and glucose-responsive transcription factors that regulate five or more metabolic enzymes. Third, we reduced the interlayer regulatory connections from the Metabolite layer to the Pathway layer by removing the interlayer regulatory connections that regulated fewer than five metabolic reactions.
Implementation. Statistical tests, clustering analysis, enrichment analysis, and transomic network analysis were done using MATLAB 2020a (The Mathworks Inc.). Visualization of transomic network in the Graph Modeling Language formats was done using Python 2.7 and VANTED 58 .

Data availability
Sequencing data measured in this study have been deposited in the DNA Data Bank of Japan Sequence Read Archive (DRA) (www. ddbj. nig. ac. jp/) under the accession no. DRA008416. All other data needed to evaluate the conclusions in the paper are present in the paper or Supplementary Materials. The code used for the analysis in this paper is available upon request.